Project Overview

Extreme storms can rapidly devastate coral reefs, yet their consequences for apex predators remain unclear. We combined year-long acoustic monitoring of 25 tagged predators (giant trevally Caranx ignobilis, Galapagos sharks Carcharhinus galapagensis, tiger sharks Galeocerdo cuvier) at French Frigate Shoals (Hawaiʻi) with bulk and compound-specific amino-acid stable isotope analyses of pre-storm muscle tissue to assess survival, space use, and trophic pathways surrounding Category-4 Hurricane Walaka (October 2018).

Image. Giant trevally (Caranx ignobilis) in the Northwestern Hawaiian Islands.
Image. Giant trevally (Caranx ignobilis) in the Northwestern Hawaiian Islands.

Data

Data here consists of bulk and amino acid δ13C and δ15N isotope values for Giant trevally (Caranx ignobilis) and Galapagos sharks (Carcharhinus galapagensis) from the Northwestern Hawaiian Islands, at French Frigate Shoals atoll (’ōlelo Hawai’i: Lalo or Mokupāpapa). We also include data used to calculate residence index using telemetry, as well as info on tagging dates and consumer fork length. Lastly, we incorporate the data from Skinner et al. (2021) reporting δ13C isotope values for essential amino acids (EAAs) collected inshore and offshore of coral reef fishes and invertebrates (plankton and squid) in the Maldives. These are used as trianing data to assign fish consumers and determine the proportion support of foodweb support to consumers.

Skinner, C., Mill, A. C., Fox, M. D., Newman, S. P., Zhu, Y., Kuhl, A., & Polunin, N. V. C. (2021). Offshore pelagic subsidies dominate carbon inputs to coral reef predators. In Sci. Adv (Vol. 7).

Linear Discriminant Analysis

Subset the data and use ‘sources’ to generate training data that will later be used to assign consumer. It is important that sources are distinct and have high accuracy in their assignment, otherwise they may not be informative as distinct energy channels/prey items.

Source groups consisted of coral reef consumers, including fishes, plankton, and squid, which reflect the primary producers of their respective habitats (benthic versus pelagic) and food webs (inshore versus offshore). Three source groups were used: “reef-associated plankton” represented by diurnal (Caesio varilineata, Caesio xanthonota) and nocturnal fishes (Myripristis violacea); “pelagic plankton” derived from squid (Uroteuthis duvaucelii); and “coral” derived from coralivorous fish (Chaetodo meyeri) found in Skinner et al. (2021).

########## LDA setup

raw.df$SampleID<-as.factor(raw.df$SampleID)
raw.df$Group<-as.factor(raw.df$Group)

# remove all samples, keep source
prod <- raw.df[!(raw.df$Group == "GT" |raw.df$Group == "GS"),]
prod$Group<-droplevels(prod$Group) # drop sample levels

#Run an LDA with a jackknifing model fit to look at error rate
#'CV = TRUE' makes the LDA run a jackknifed - leave one out - model fit
All.lda <- lda(Group ~Thr + Val + Leu + Phe + Lys, data = prod, CV = TRUE)

Here, show LDA classification of sources performed with 100% success.

#Create a table which compares the classification from the LDA model to the actual classification
All.reclass <- table(prod$Group, All.lda$class) #100% success

#Total percent of samples correctly classified is the sum of the diagonal of this table
sum.dig<-sum(diag(prop.table(All.reclass))) #100% !

accur<-mean(All.lda$class == prod$Group) # mean accuracy

#Percent of each producer Group correctly reclassified: 100% for all groups
reclass<-diag(prop.table(All.reclass, 1))

# print results
All.reclass
##                  
##                   Coral PelagicPlankton ReefPlankton
##   Coral               9               0            0
##   PelagicPlankton     0               8            0
##   ReefPlankton        0               0           25

Show the training data with the coefficients of linear discriminants and proportion of trace for LD1 and LD2. This shows the EAA that are driving patterns across the data visualization, as well as the % trace explained on each axis. Here 79% of variation is explained by PC1 and this is strongly affected by threonine (Thr) and phenylalanine (Phe).

#Create a training LDA function from the library data
All.train <- lda(Group~ Thr + Val + Leu + Phe + Lys, data = prod)
All.train
## Call:
## lda(Group ~ Thr + Val + Leu + Phe + Lys, data = prod)
## 
## Prior probabilities of groups:
##           Coral PelagicPlankton    ReefPlankton 
##       0.2142857       0.1904762       0.5952381 
## 
## Group means:
##                        Thr       Val       Leu       Phe       Lys
## Coral            -8.378889 -21.12333 -19.74556 -15.02444 -12.83556
## PelagicPlankton  -8.485000 -22.72750 -23.35125 -24.77125 -16.25500
## ReefPlankton    -13.307200 -24.07960 -23.55920 -21.27280 -18.21480
## 
## Coefficients of linear discriminants:
##              LD1         LD2
## Thr  1.226787772 -0.47940410
## Val -0.071830357  0.66667233
## Leu -0.269481800 -0.14587473
## Phe -1.480860015  0.09941465
## Lys  0.009195137 -1.08451168
## 
## Proportion of trace:
##    LD1    LD2 
## 0.7913 0.2087
# save as txt file
sink("output/Walak_LDA_Train_Info.txt")
print(All.train)
sink()

Classify consumers (giant trevally or Galapagos sharks) as belonging to the three source groups.

  • We show the classification gives about 1/2 of the trevally as inshore vs. offshore, with 3 of 4 sharks being inshore.
####### ######## ######## ######## #####
######## Classify Trevaly and Sharks ####
######## ######## ######## ######## #####

# remove all samples, keep source
cons.df <- raw.df[(raw.df$Group == "GT" |raw.df$Group == "GS"),]
cons.df$Group<-droplevels(cons.df$Group) # drop sample levels

cons.predict <- predict(object = All.train, newdata = cons.df)

cons.predict.data <- data.frame(SampleID = cons.df$SampleID, Group = cons.df$Group, LD.class = cons.predict$class, cons.predict$x)

All.reclass <- table(cons.df$Group, cons.predict$class) 
All.reclass
##     
##      Coral PelagicPlankton ReefPlankton
##   GS     0               1            3
##   GT     0               7            6
#  GS: 1 as PelagicPlank, 3 as ReefPlankton
#  GT: 7 as PelagicPlank, 6 as ReefPlankton

Plot the LDA with producer groups and the consumers.

#### plot it
raw.norm.LD.data$Group<-factor(raw.norm.LD.data$Group, 
                      levels=c("PelagicPlankton", "ReefPlankton", "Coral", 
                               "GT", "GS"))

lda.5.colors <- c("dodgerblue", "springgreen4","coral", "gold", "darkgoldenrod")


LDA.scat.Group <- ggplot(raw.norm.LD.data, aes(x = LD1, y = LD2))+
  geom_point(size = 3.5, color="black",
             aes(fill = Group, color = Group, shape = Group), alpha = 0.9) +
  stat_ellipse(type = "norm", level = 0.90, linewidth = 0.5, aes(lty = Group, color = Group))+
  scale_color_manual(values = lda.5.colors)+
  scale_fill_manual(values = lda.5.colors)+
  scale_shape_manual(values = c(24, 24, 24, 21, 21))+
  scale_linetype_manual(values = c(2,2,2,0,0))+
  xlab("LD1 (79.1%)") + ylab("LD2 (20.9%)") +
  ggtitle(expression(paste("LDA: Measured-", delta^{13}, C[EAA]))) +
  theme_classic() + guides(lty = "none")

LDA.scat.Group
**Figure 6.** Linear discriminate analysis (LDA) of δ^13^C~EAA~ fingerprints from reef-associated plankton, pelagic plankton and coral compared with those of predators. Ellipses denote 95% confidence intervals for each source group.

Figure 6. Linear discriminate analysis (LDA) of δ13CEAA fingerprints from reef-associated plankton, pelagic plankton and coral compared with those of predators. Ellipses denote 95% confidence intervals for each source group. 

dev.copy(pdf, "figures/Walak_LDA_sorc_consumer.pdf", width=7, height=6)
dev.off()

MixSIAR

LDA is great at classification, but cannot handle mixtures. Use a 3 member mixing model (reef plankton, pelagic plankton, coral) in MixSIAR to generate estimates of resource use. Once these models run, then perform a posteriori aggregation to generate estimates for 2 sources: inshore vs. offshore.

First, load in sources and generate the means and SD for each normalized-EAA.

################################################
## MIXSIAR #####################################
################################################

# create a new data frame with desired columns for the consumers (mixtures)
Mix.df.EAA <- raw.norm.LD.data %>%
  dplyr::select("SampleID", "Group", "Thrn", "Valn", "Leun", "Phen", "Lysn")

# remove sources
Mix.df.EAA.consm <- Mix.df.EAA[!(Mix.df.EAA$Group == "ReefPlankton" | Mix.df.EAA$Group == "PelagicPlankton" |
                                  Mix.df.EAA$Group == "Coral"),]
Mix.df.EAA.consm$Group <- droplevels(Mix.df.EAA.consm$Group)
Mix.df.EAA.consm<-na.omit(Mix.df.EAA.consm)

# write and export a csv of the data frame
write.csv(Mix.df.EAA.consm, "output/EAAn_Mix/Walak_Mix.df.EAAn.consumer.csv", row.names = FALSE)

###########
# SOURCES #
###########

# NORMALIZED DATA
# Reef and pelagic plankton, coral

# create a new data frame with desired columns for the primary producers (sources)
Source.df.EAA <-  Mix.df.EAA[(Mix.df.EAA$Group == "ReefPlankton" | Mix.df.EAA$Group == "PelagicPlankton" |
                               Mix.df.EAA$Group == "Coral"),]

Source.df.EAA <- Source.df.EAA %>%
  dplyr::select("Group", "Thrn", "Valn", "Leun", "Phen", "Lysn")

##### source mean and SD
### Thr mean and SD
# calculate mean Thr values for the sources
mean.Thr.df.EAA <- aggregate(Thrn ~ Group, data = Source.df.EAA, FUN = mean)
colnames(mean.Thr.df.EAA) <- c("Group", "Mean.Thrn")

# calculate standard deviations of Thr values for the sources
sd.Thr.df.EAA <- aggregate(Thrn ~ Group, data = Source.df.EAA, FUN = sd)
colnames(sd.Thr.df.EAA) <- c("Group", "SD.Thrn")

### Val mean and SD
# calculate mean Val values for the sources
mean.Val.df.EAA <- aggregate(Valn ~ Group, data = Source.df.EAA, FUN = mean)
colnames(mean.Val.df.EAA) <- c("Group", "Mean.Valn")

# calculate standard deviations of Val values for the sources
sd.Val.df.EAA <- aggregate(Valn ~ Group, data = Source.df.EAA, FUN = sd)
colnames(sd.Val.df.EAA) <- c("Group", "SD.Valn")

### Leu mean and SD
# calculate mean Leu values for the sources
mean.Leu.df.EAA <- aggregate(Leun ~ Group, data = Source.df.EAA, FUN = mean)
colnames(mean.Leu.df.EAA) <- c("Group", "Mean.Leun")

# calculate standard deviations of Leu values for the sources
sd.Leu.df.EAA <- aggregate(Leun ~ Group, data = Source.df.EAA, FUN = sd)
colnames(sd.Leu.df.EAA) <- c("Group", "SD.Leun")

### Phe mean and SD
# calculate mean Phe values for the sources
mean.Phe.df.EAA <- aggregate(Phen ~ Group, data = Source.df.EAA, FUN = mean)
colnames(mean.Phe.df.EAA) <- c("Group", "Mean.Phen")

# calculate standard deviations of Phe values for the sources
sd.Phe.df.EAA <- aggregate(Phen ~ Group, data = Source.df.EAA, FUN = sd)
colnames(sd.Phe.df.EAA) <- c("Group", "SD.Phen")

### Lys mean and SD
# calculate mean Phe values for the sources
mean.Lys.df.EAA <- aggregate(Lysn ~ Group, data = Source.df.EAA, FUN = mean)
colnames(mean.Lys.df.EAA) <- c("Group", "Mean.Lysn")

# calculate standard deviations of Phe values for the sources
sd.Lys.df.EAA <- aggregate(Lysn ~ Group, data = Source.df.EAA, FUN = sd)
colnames(sd.Lys.df.EAA) <- c("Group", "SD.Lysn")

#### sample size
# calculate the number of samples for each of the sources
df.size.EAA <- aggregate(Lysn ~ Group, data = Source.df.EAA, FUN = length)
colnames(df.size.EAA) <- c("Group", "Lys.size")

# bind the values calculated above into one data frame
source.agg.df.EAA <- cbind(df.size.EAA, 
                           mean.Thr.df.EAA[2], sd.Thr.df.EAA[2], 
                           mean.Val.df.EAA[2], sd.Val.df.EAA[2], 
                           mean.Leu.df.EAA[2], sd.Leu.df.EAA[2], 
                           mean.Phe.df.EAA[2], sd.Phe.df.EAA[2], 
                           mean.Lys.df.EAA[2], sd.Lys.df.EAA[2])
colnames(source.agg.df.EAA) <- c("", "n", "MeanThrn", "SDThrn", "MeanValn", "SDValn", 
                                 "MeanLeun", "SDLeun", "MeanPhen", "SDPhen", "MeanLysn", "SDLysn")

# write and export a csv of the source data frame
write.csv(source.agg.df.EAA, "output/EAAn_Mix/Walak_source.agg.df.EAAn.csv", row.names = FALSE, na = "")

Now, load in the consumers. We use a TDF of 0 (with SD of 0) assuming direct routing of EAAs to consumers.

# load consumer data and assign factors
# SampleID is a fixed factor
cons.EAAn <- load_mix_data(filename = "output/EAAn_Mix/Walak_Mix.df.EAAn.consumer.csv",
                          iso_names = c("Thrn", "Valn", "Leun", "Phen", "Lysn"),
                          factors = "SampleID",
                          fac_random = FALSE,
                          fac_nested = FALSE,
                          cont_effects = NULL)

# load source data
source.EAAn <- load_source_data(filename = "output/EAAn_Mix/Walak_source.agg.df.EAAn.csv",
                               source_factors = NULL,
                               conc_dep = FALSE,
                               data_type = "means",
                               cons.EAAn)

#Load TDF data
discr.EAAn <- load_discr_data(filename = "output/EAAn_Mix/MixSIAR_EAAn_TDFs.csv", cons.EAAn)

Set up here for the model (JAGs), run the model and save the output as an R.data file (can be memory intensive). For the sake of the Rmd-knitr, this is set to eval=FALSE.

################## ################## ################# 
################## RUN MIXSIAR MODEL BY SAMPLEID #######

# define the model and error structure and write JAGS model file
# PROCESS ERROR ONLY - HAVE TO DO THIS WHEN RUNNING MIXTURE POINTS INDIVIDUALLY

model_filename.EAA <- "output/MixSIAR_models/EAAn_SampleID_model/MixSIAR_Model_SampleID_fixed.txt"
resid_err <- FALSE
process_err <- TRUE
write_JAGS_model(model_filename.EAA, resid_err, process_err, cons.EAAn, source.EAAn)

################################################
############# run the JAGS model

set.seed(12)
consm.jags.mod.EAAn <- run_model(run = "normal", cons.EAAn, source.EAAn, discr.EAAn, model_filename.EAA, alpha.prior = 1, resid_err, process_err)

#save RDS object of JAGS model
saveRDS(consm.jags.mod.EAAn, file = "output/MixSIAR_models/EAAn_SampleID_model/MixSIAR_Model_SampleID_rjags_norm.rds")
# code below will read back in the RDS object
consm.jags.mod.EAAn <- readRDS(file = "output/MixSIAR_models/EAAn_SampleID_model/MixSIAR_Model_SampleID_rjags_norm.rds")

Generate plots for MixSIAR diagnostics. See the MixSIAR model statistics here. Gelman-Rubin values are < 1.05, and 2 of 3 chains in Geweke diagnostic perform well at <5% z-score.

# note the processing here opens files in quartz and results in an issue with knitting. This is why code chunk is set to "eval = FALSE". Code runs otherwise.

# process diagnostics, summary stats, and posterior plots
output_options_EAAn <- list(
  summary_save = TRUE,                   
  summary_name = "output/MixSIAR_models/EAAn_SampleID_model/MixSIAR_SampleID_summary_stats.EAAn",
  sup_post = FALSE,                       
  plot_post_save_pdf = FALSE,              
  plot_post_name = "output/MixSIAR_models/EAAn_SampleID_model/MixSIAR_SampleID_posterior_density.EAAn",   
  sup_pairs = TRUE,                       
  plot_pairs_save_pdf = FALSE,             
  plot_pairs_name = "output/MixSIAR_models/EAAn_SampleID_model/MixSIAR_SampleID_pairs_plot.EAAn",         
  sup_xy = FALSE,                          
  plot_xy_save_pdf = TRUE,                
  plot_xy_name = "output/MixSIAR_models/EAAn_SampleID_model/MixSIAR_SampleID_xy_plot.EAAn",             
  gelman = TRUE,                         
  heidel = FALSE,                         
  geweke = TRUE,                         
  diag_save = TRUE,                       
  diag_name = "output/MixSIAR_models/EAAn_SampleID_model/MixSIAR_SampleID_diagnostics.EAAn",
  indiv_effect = FALSE,
  plot_post_save_png = FALSE,            
  plot_pairs_save_png = FALSE,            
  plot_xy_save_png = FALSE,
  graphics.off())
output_diagnostics(consm.jags.mod.EAAn, cons.EAAn, source.EAAn, output_options_EAAn)
## 
## ################################################################################
## # Gelman-Rubin Diagnostic
## ################################################################################
## 
## Generally the Gelman diagnostic should be < 1.05
## 
## Out of 106 variables: 10 > 1.01
##                       0 > 1.05
##                       0 > 1.1
## 
## The worst variables are:
##              Point est. Upper C.I.
## p.fac1[17,1]   1.016869   1.035497
## p.fac1[9,1]    1.015209   1.039946
## p.fac1[7,1]    1.014457   1.032430
## loglik[9]      1.014359   1.024013
## loglik[16]     1.012579   1.022181
## p.fac1[16,1]   1.012560   1.027758
## p.fac1[1,1]    1.011651   1.037135
## p.global[1]    1.011651   1.037135
## p.fac1[12,3]   1.010392   1.039600
## p.fac1[4,1]    1.010207   1.030760
## 
## ################################################################################
## # Geweke Diagnostic
## ################################################################################
## 
## The Geweke diagnostic is a standard z-score, so we'd expect 5% to be outside +/-1.96
## Number of variables outside +/-1.96 in each chain (out of 106):
## 
##        Chain 1 Chain 2 Chain 3
## Geweke      23       3       1

Output for figures causes quartz device to open, which also gives issue with knit in generating the html output. For the sake of the Rmd-knitr, this is set to eval=FALSE.

# set max to get all data
options(max.print = 3000)
out.JAGS <- output_JAGS(consm.jags.mod.EAAn, cons.EAAn, source.EAAn, output_options_EAAn)

# get the output stats for the model
# output_stats(consm.jags.mod.EAAn, cons.EAAn, source.EAAn, output_options_EAAn)

Perform source aggregation to “inshore vs. offshore”.

######### a posteriori aggregation of sources
# combine benthic algae and POM into one autochthonous source
consm.in.out.combined.EAAn <- combine_sources(consm.jags.mod.EAAn, cons.EAAn, source.EAAn, alpha.prior = 1,
                                          groups = list(Inshore = c("ReefPlankton", "Coral"),
                                                        Offshore = "PelagicPlankton"))

# print and save the interval plot
plot_intervals(consm.in.out.combined.EAAn, toplot = "fac1")
dev.copy(pdf, "output/MixSIAR_models/EAAn_SampleID_model/combines_output_intervals.EAAn.pdf", width=10, height=10)
dev.off()

#### # summary statistics for a posteriori aggregation of sources (inshore vs offshore)
# ***** NOTE ***** you will need to create a blank txt file if you run and re-run the aggregation. The file is not overwritten when you re-save in summary_stat(...). It will add outputs sequentially, which will make life terrible and code not-reproducible.

# make blank txt if re-running summary_stat(...)
file.create("output/MixSIAR_models/EAAn_SampleID_model/MixSIAR_SampleID_summary_stats_a_posteriori.EAAn.txt")

# run summary_stat and save output
summary_stat(consm.in.out.combined.EAAn, toprint = "fac1", meanSD = TRUE, savetxt = TRUE,
             quantiles = c(0.025, 0.05, 0.25, 0.5, 0.75, 0.95, 0.975),
             filename = "output/MixSIAR_models/EAAn_SampleID_model/MixSIAR_SampleID_summary_stats_a_posteriori.EAAn")
################################################
### CLEAN MIXSIAR MODEL BY SAMPLEID RESULTS ###
################################################

SAMP.mixSIAR.EAAn.output.3 <- read.table(
  "output/MixSIAR_models/EAAn_SampleID_model/MixSIAR_SampleID_summary_stats.EAAn.txt",
  sep = '\t', header = FALSE, skip = 6)

SAMP.mixSIAR.EAAn.output.3 <- as.data.frame(SAMP.mixSIAR.EAAn.output.3)
colnames(SAMP.mixSIAR.EAAn.output.3) <- "full.data"

# remove all extra spaces using "squish", make them 1 space only
library(stringr)
SAMP.mixSIAR.EAAn.output.3$full.data <- str_squish(SAMP.mixSIAR.EAAn.output.3$full.data)

# rename columns
library(plyr)
samp.mix.EAAn.out.3 <- separate(SAMP.mixSIAR.EAAn.output.3, full.data,  
                               into = c("factor", "Mean", "SD", "2.5.perc", "5.perc", "25.perc", "50.perc",
                                        "75.perc", "95.perc", "97.5.perc"), sep=" ")

# parse the factors here into 3 columns
new.cols.EAAn.3 <- stringr::str_split_fixed(samp.mix.EAAn.out.3$factor, "\\.", 3) %>%
  as.data.frame() %>%
  setNames(c("proportion", "SampleID", "Source"))

new.cols.EAAn.3 <- new.cols.EAAn.3 %>%
  dplyr::select("SampleID", "Source")

# combine the new, separated columns with the data, verify the levels match, then remove "factor"
Samp.out.cleaned.EAAn.3 <- cbind(new.cols.EAAn.3, samp.mix.EAAn.out.3[1:10])
Samp.out.cleaned.EAAn.3 <- Samp.out.cleaned.EAAn.3 %>%
  dplyr::select(-factor)

Samp.out.cleaned.EAAn.3<-Samp.out.cleaned.EAAn.3[-1,]

write.csv(Samp.out.cleaned.EAAn.3, "output/EAAn_Mix/Walak_MixSIAR_SampleID.df.EAAn.3sources.csv", row.names = FALSE)

3 source MixSIAR plot

Plot of three sources and their contribution to consumers.

################## ready to plot ################## 

# bring in some metadata to make sense of it
Mix.samp.EAAn.3 <- merge.data.frame(raw.norm.LD.data, Samp.out.cleaned.EAAn.3, by="SampleID")
Mix.samp.EAAn.3$Source <- factor(Mix.samp.EAAn.3$Source, 
                                 levels = c("PelagicPlankton", "ReefPlankton", "Coral"))
Mix.samp.EAAn.3$Mean<-as.numeric(Mix.samp.EAAn.3$Mean) 

# plot it: some general formatting here for the plot...
Fig.formatting <- 
  theme(axis.ticks.length = unit(0.2, "cm"),
        axis.text.x = element_text(size = 12, color = "black"),
        axis.text.y = element_text(size = 12, color = "black"),
        axis.title.x = element_text(size = 14, color = "black"),
        axis.title.y = element_text(size = 14, color = "black"),
        legend.text = element_text(size = 10, color = "black"),
        legend.title = element_blank(),
        axis.ticks = element_line(color = "black"),
        legend.key.size = unit(0.6, "cm"))


## plot the 3 sources and consumers
Mix.samp.plot.EAAn.3 <- ggplot(data = Mix.samp.EAAn.3, aes(x = Group, y = Mean, fill = Source)) +
  geom_boxplot(alpha = 0.9, color = "black")+
  geom_point(aes(fill = Source), pch = 21, alpha = 0.5, color = "black", 
             position = position_dodge(0.75))+
  scale_fill_manual(values = c("dodgerblue", "springgreen4","coral")) +
  ylab("Proportional Contribution") +
  xlab(NULL) +
  ylim(0,1) +
  guides(fill = guide_legend(title = expression(paste(AA[ESS], " ", Source)))) +
  theme_classic() + Fig.formatting + guides(color = "none")

#print and save
Mix.samp.plot.EAAn.3
**Figure 7a.** Results of the MixSIAR Bayesian mixing model using δ^13^C~EAA~ fingerprints. Proportional contribution of pelagic plankton, reef-associated plankton and corals to Giant trevally and Galapagos sharks.

Figure 7a. Results of the MixSIAR Bayesian mixing model using δ13CEAA fingerprints. Proportional contribution of pelagic plankton, reef-associated plankton and corals to Giant trevally and Galapagos sharks. 

dev.copy(pdf, "figures/Walak_Fig.mixsiar.3sour.EAAn.pdf", width = 6, height = 7)
dev.off()

2 member MixSIAR plot

Using aggregated sources, plot the % contribution of inshore vs. offshore food web support to consumers

samp.mix.df<-Samp.out.cleaned.EAAn
samp.mix.reduced <- samp.mix.df %>% 
  dplyr::select(SampleID, Source, Mean, SD)

################# ################### ################ ################

# bring in some metadata to make sense of it
Mix.samp.EAAn.2 <- merge.data.frame(raw.norm.LD.data, samp.mix.reduced, by="SampleID")
Mix.samp.EAAn.2$Source <- factor(Mix.samp.EAAn.2$Source, 
                                 levels = c("Offshore", "Inshore"))
Mix.samp.EAAn.2$Mean<-as.numeric(Mix.samp.EAAn.2$Mean) 

write.csv(Mix.samp.EAAn.2, "output/Walak_Mix.samp.EAAn.2sources.csv")

########### plot
Mix.samp.plot.EAAn.InOff <- ggplot(data = Mix.samp.EAAn.2, aes(x = Group, y = Mean, fill = Source)) +
  geom_boxplot(alpha = 0.9, color = "black")+
  geom_point(aes(fill = Source), pch = 21, alpha = 0.5, color = "black", position = position_dodge(0.75))+
  scale_fill_manual(values = c("dodgerblue","springgreen4")) +
  ylab("Proportional Contribution") +
  xlab(NULL) +
  ylim(0,1) +
  guides(fill = guide_legend(title = expression(paste(AA[ESS], " ", Source)))) +
  theme_classic() + Fig.formatting + guides(color = "none")

#print and save grouped with 2 categories
Mix.samp.plot.EAAn.InOff
**Figure.** Results of the MixSIAR Bayesian mixing model using δ^13^C~EAA~ fingerprints. Proportional contribution of inshore versus offshore resources to giant trevally and Galapagos sharks.

Figure. Results of the MixSIAR Bayesian mixing model using δ13CEAA fingerprints. Proportional contribution of inshore versus offshore resources to giant trevally and Galapagos sharks. 

dev.copy(pdf, "figures/Walak_Fig.mixsiar.2sour.EAAn.pdf", width = 4, height = 7)
dev.off()

Residency index

Her we use the residency index, or the events where the consumers were recorded by receiver arrays at Fish Frigate Shoals relative to the days the receivers were recording, to compare the proportion of inshore food web support to consumers.

We show that there is a significant correlation between the residency index and the % inshore food web use, dominated by the Giant trevlly, which are best represented in our dataset. The relationship is significant R2 = 0.55 and p<0.001.

Inshore.df<-Mix.samp.EAAn.2[(Mix.samp.EAAn.2$Source=="Inshore"),]
# write.csv(Inshore.df, "output/Inshore.df.csv")

# group mean model
Inshore.mod<-lm(Mean~Group, data=Inshore.df) #p=0.610, no difference in % inshore by group
Group.model<-anova(Inshore.mod)

# residency by mean model
Res.Inshore.mod<-lm(Mean~Residency.index, data=Inshore.df)
Res.sum<-summary(Res.Inshore.mod) #R2-adj = 0.55,  #p<0.001
anova(Res.Inshore.mod)
## Analysis of Variance Table
## 
## Response: Mean
##                 Df   Sum Sq Mean Sq F value    Pr(>F)    
## Residency.index  1 0.096760 0.09676  20.943 0.0003636 ***
## Residuals       15 0.069302 0.00462                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This relationship is unchanged by the presence or absence of sharks in the data. Here, just giant trevally are used in the model, Adj-R2 = 0.56 and p=0.002.

# residency model if GS removed
Res.Inshore.mod.GT<-lm(Mean~Residency.index, data=Inshore.df[(Inshore.df$Group=="GT"),])
Res.sum<-summary(Res.Inshore.mod.GT) #R2-adj = 0.56, p=0.002
anova(Res.Inshore.mod.GT)
## Analysis of Variance Table
## 
## Response: Mean
##                 Df  Sum Sq  Mean Sq F value  Pr(>F)   
## Residency.index  1 0.08742 0.087420  16.343 0.00194 **
## Residuals       11 0.05884 0.005349                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residency by proportion of inshore foodweb support

Plot the relationship for residency index for all sames (sharks and trevally).

Inshore.prop.plant<-ggplot(Inshore.df, aes(x=Residency.index, y=Mean)) +
  geom_smooth(method = lm, color="seagreen", alpha=0.2, fill="mediumseagreen")+
  geom_point(color="seagreen", size=2.5, aes(shape=Group))+
  ylab("Proportion Inshore-EAA")+
  xlab("Residency Index")+
  annotate("text", x = 0.08, y = 0.68, label = "Adj-R2=0.55, p<0.001", size = 3) +
  theme_classic()

Inshore.prop.plant
**Figure 7B.** Relationship between the proportion of inshore carbon sources utilized versus and the predator residency index. The correlation is significant (ANOVA, f = 20.943, p <0.001)

Figure 7B. Relationship between the proportion of inshore carbon sources utilized versus and the predator residency index. The correlation is significant (ANOVA, f = 20.943, p <0.001) 

ggsave("figures/Walak_Inshore.prop.EAAn.pdf", width=5, height=7)